\chapter{Dust Formation}
\label{chap:dust}

In Chapter \ref{chap:chemical} we discussed the chemical forms into which the
outflows would condensate if the system becomes an equilibrium. But for the
high speed of the outflows it is difficult to get into real equilibrium. So
here we study the dynamic dust formation, assuming single type of atom for
simplicity.

%===============================================================================
\section{Abundance Change}
%===============================================================================

For a reaction:
\[
1 + 2 \longleftrightarrow 3 + \gamma
\]
we can write the number density change of species 1 as:
\begin{equation}
\frac{{\rm d}n_1}{{\rm d}t} = 
  -\lambda_{12} \, n_1 n_2 + \lambda_\gamma n_3
\label{eq:number_density}
\end{equation}
where $\lambda_{12}$ is the reaction rate in the unit of $cm^3s^{-1}$ and 
$\lambda_\gamma$ is the photo-dissociation rate of species 3 with the unit of
$s^{-1}$.

In a expanding or a collapsing system the number densities change with volume. 
We'd like to study a quantity that does not depend on volume directly so that
we use the fraction of each species in number instead of number density:
\[
Y_i \equiv n_i / n_{total}
\]
where $Y_i$ is called the abundance of species i, $n_i$ is the number density
of species i and $n_{total}$ is the number density of the unit particles, here
the number density of total atoms. So $Y_i$ here indicates the number
of species i per atom. And
\[
n_{total} = n_{atom} = n_{baryon} \times Y_{atom} 
\]
where $n_{baryon} = \rho N_A$ is the baryon number density and $Y_{atom}$ is
the atom fraction, i.e. number of atoms per baryon. Take pure $^{12}$C for 
example, $Y_{atom} = 1 / 12$. So Eq.\ref{eq:number_density} can be rewrite as
\begin{equation}
\frac{{\rm d}Y_1}{{\rm d}t} = 
  -\lambda_{12} n_{atom} Y_1 Y_2 + \lambda_\gamma Y_3 =
  -\lambda_{12} \rho N_A Y_{atom} Y_1 Y_2 + \lambda_\gamma Y_3
\label{eq:abundance}
\end{equation}
where $\rho$ is the mass density and $N_A$ is the Avogadro number. Then what
we need to study is how the abundances of each species change with time.

%===============================================================================
\section{Condensation Effective Rates}
%===============================================================================

For a single kind of atom, we study the effective condensation rate here.
The condensation here means one atom attaches to an existing multi-atom
particle:
\begin{equation}
P_1 + P_i \to P_{i+1} + \gamma
\label{eq:reaction}
\end{equation} 
where $P_i$ denotes a particle with $i$ atoms in it. Here we are neglecting
the photo-dissociation for large stable particles.

Let's assume that the condensation rate is proportional to the surface area
(or cross section) of the larger particle, which goes like $N^{2/3}$, i.e.
\[
\lambda \propto A \propto N^{2/3}
\]
or
\[
\lambda_i = \lambda_1 N_i^{2/3}
\]
where $\lambda$ is the rate, $A$ is the surface area and $N$ is the number 
of atoms in the particle.

To make the calculation of particles with a large atom number possible, we 
introduce
the ``bin" after a continuous network. There will a number of particles in a 
bin and we assume \textbf{steady state} within a bin.
The rate out of a bin is like
\[
\lambda_{bin,last}Y_{bin,last} = \lambda^{eff}_{bin}Y_{bin}
\]
\[
Y_{bin} = \sum_{i\in bin} Y_i
\]
\[
\lambda^{eff}_{bin} = \lambda_{bin,last} \left(
  \frac{Y_{bin,last}}{Y_{bin}}
  \right)_{ss}
  = \lambda_{bin,last} \left(
  \frac{1}{\sum_{i\in bin} \frac{Y_i}{Y_{bin,last}} }
  \right)_{ss}
\]
For steady state,
\[
\frac{Y_i}{Y_j} = \frac{\lambda_i}{\lambda_j}
\]
so,
\[
\lambda^{eff}_{bin} = \lambda_{bin,last} \left(
  \frac{1}{\sum_{i\in bin} \frac{\lambda_{bin,last}}{\lambda_i} }
  \right)
  = \frac{1}{\sum_{i\in bin} \frac{1}{\lambda_i} }
  = \lambda_1 \frac{1}{\sum_{i\in bin} N_i^{-2/3}}
\]
If we apply the integral approximation,
\begin{equation}
\lambda^{eff}_{bin}
  \approx \lambda_1 \left( \int_{bin,fist}^{bin,last} x^{-2/3} {\rm d} x
  \right)^{-1}
  = \lambda_1 \frac{1}{N_{bin,last}^{1/3} - N_{bin,first}^{1/3}}
\label{eq:rate}
\end{equation}

We estimate $\lambda_1$ as two single atoms collide together with thermal 
velocity and the transverse section.
\begin{equation}
\lambda_1 = \: <\! \sigma_1 v\! > \: = \pi R^2 \sqrt{3k_BT/m}
\label{eq:k1}
\end{equation}
where $\sigma_1$ is the cross section for a single atom and $R$ is the radius
of a single atom, $v$ is the thermal velocity which is taken to be the average
value in 3 dimensions, $k_B$ is the Boltzmann constant, $T$ is the temperature
and $m$ is the mass of a single atom.

%===============================================================================
\section{Atom Conservation}
%===============================================================================

The abundance change rates of $Y_{bin,i}$, $Y_{bin,i+1}$ and $Y_1$
due to 
reaction ${\rm Bin}_i + X P_1 \to {\rm Bin}_{i+1} + \gamma$ are
\[
\frac{{\rm d}Y_{bin,i}}{{\rm d}t} =
  -\lambda_i^{eff} n_{atom} Y_{bin,i} Y_1
\]
\[
\frac{{\rm d}Y_{bin,i+1}}{{\rm d}t} =
  \lambda_i^{eff} n_{atom} Y_{bin,i} Y_1
\]
\[
\frac{{\rm d}Y_1}{{\rm d}t} =
  X( -\lambda_i^{eff} n_{atom} Y_{bin,i} Y_1 )
\]
where $X$ is the factor to be determined that ensure
the total atom number conserved. Atom conservation requires
\[
N_{bin,i} \times \frac{{\rm d}Y_{bin,i}}{{\rm d}t} +
N_{bin,i+1} \times \frac{{\rm d}Y_{bin,i+1}}{{\rm d}t} +
\frac{{\rm d}Y_1}{{\rm d}t} = 0.
\] 
So we have
\[
X = N_{bin,i+1} - N_{bin,i}
\]

Similarly the abundance change rates of $Y_{last}$, $Y_{bin,1}$ and 
$Y_1$ due to reaction $P_{network,last} + X P_1 \to {\rm Bin}_1 + \gamma$ are
\[
\frac{{\rm d}Y_{last}}{{\rm d}t} =
  -\lambda_{last} n_{atom} Y_{last} Y_1
\]
\[
\frac{{\rm d}Y_{bin,1}}{{\rm d}t} =
  \lambda_{last} n_{atom} Y_{last} Y_1
\]
\[
\frac{{\rm d}Y_1}{{\rm d}t} =
  (-\lambda_{last} n_{atom} Y_{last} Y_1) \times
  (N_{bin,1} - N_{last})
\]

\section{Network Matrix}
Let's look at the reaction between the last species in the continuous network
and the first bin, 
the changes of the root-finding functions are (how to call them?)
(denoting $n_{atom}$ as $n_a$):
\[
f_{last} \: +\!\!= \: \lambda_{last} \, n_a Y_{last} Y_1
\]
\[
f_{1} \: +\!\!= \: (\lambda_{last} \, n_a Y_{last} Y_1) 
  (N_{bin,1} - N_{last})
\]
To find the abundance changes in the continuous network, we need to solve the
matrix
\[
AX=B
\]
The matrix element changes due to this reaction are
\[
A_{last,last} \: +\!\!= \: \lambda_{last} \, n_a Y_1
\]
\[
A_{last,1} \: +\!\!= \: \lambda_{last} \, n_a Y_{last}
\]
\[
A_{1,last} \: +\!\!= \: (\lambda_{last} \, n_a Y_1)(N_{bin,1} - N_{last})
\]
\[
A_{1,1} \: +\!\!= \: (\lambda_{last} \, n_a Y_{last})(N_{bin,1} - N_{last})
\]
The right-hand-side vector changes are
\[
B_{last} \: +\!\!= \: -f_{last} = -\lambda_{last} \, n_a Y_{last} Y_1 
\]
\[
B_1 \:+\!\!=\: -f_1 = -(\lambda_{last}\,n_a Y_{last} Y_1)(N_{bin,1}-N_{last})
\]

Now let's look at the reactions between bins. 
${\rm Bin}_i + (N_{bin,i+1}-N_{bin,i})P_1 \to {\rm Bin}_{i+1}$,
where $i$ goes from 1 to $N-1$ with $N$ as the number of bins.
\[
\left( \frac{{\rm d}Y_1}{{\rm d}t}\right) _i =
  ( -\lambda^{eff}_{bin,i} \, n_a Y_{bin,i} Y_1 )(N_{bin,i+1} - N_{bin,i})
\]
\[
(f_{1})_i \: +\!\!= \: 
  (\lambda^{eff}_{bin,i} \, n_a Y_{bin,i} Y_1) (N_{bin,i+1} - N_{bin,i})
\]
\[
(A_{1,1})_i \: +\!\!= \: 
  (\lambda^{eff}_{bin,i} \, n_a Y_{bin,i}) (N_{bin,i+1} - N_{bin,i})
\]
\[
(B_1)_i \:+\!\!=\: -(f_1)_i = 
  -(\lambda^{eff}_{bin,i} \, n_a Y_{bin,i} Y_1) (N_{bin,i+1} - N_{bin,i})
\]
We first solve the matrix equation for current $Y_{bin,i}$, and then compute
the abundance changes in bins using the new abundances in the network 
(as showed below).
Then feed the new bin abundances back to the matrix. Iterate until converge. 

%===============================================================================
\section{Bin Abundances Change}
%===============================================================================

Now let's look at the abundance changes in bins for each reaction discussed
above. We are going to solve these changes implicitly. For $Y_{bin,1}$,
\[
\frac {Y_{bin,1}(t+\Delta t) - Y_{b1}(t)} {\Delta t} =
  \lambda_{last} \, n_a Y_{last}(t+\Delta t) Y_1(t+\Delta t) -
  \lambda^{eff}_{bin,1} \, n_a Y_{bin,1}(t+\Delta t) Y_1(t+\Delta t)
\]
\[
Y_{bin,1}(t+\Delta t) = \frac
  {\lambda_{last} \, n_a Y_{last}(t+\Delta t) Y_1(t+\Delta t) \Delta t +
    Y_{bin,1}(t)}
  {1 + \lambda^{eff}_{bin,1} \, n_a Y_1(t+\Delta t) \Delta t}
\]
Similarly for $Y_{bin,i}$ as $i$ goes from 2 to $N-1$:
\[
Y_{bin,i}(t+\Delta t) = \frac
  {\lambda^{eff}_{bin,i-1} \, n_a Y_{bin,i-1}(t+\Delta t) Y_1(t+\Delta t) 
    \Delta t +
    Y_{bin,i}(t)}
  {1 + \lambda^{eff}_{bin,i} \, n_a Y_1(t+\Delta t) \Delta t}
\]
For $Y_{bin,N}$ we don't have any destructive terms for now,
\[
Y_{bin,N}(t+\Delta t) = 
  \lambda^{eff}_{bin,N-1} \, n_a Y_{bin,N-1}(t+\Delta t) Y_1(t+\Delta t) 
    \Delta t +
    Y_{bin,N}(t)
\]

\begin{figure}[h!]
\centering
\includegraphics[width=\figuresize\textwidth]{figures/dust.pdf}
\caption{
Dust size distribution for a Ia trajectory.
}
\label{fig:dust}
\end{figure}

